Author

Diego Pinto

Set Up

Code
# Set working directory path
setwd("~/Desktop/R_projects/msc")

# Clear Global Environment
rm(list = ls()) 

# Load packages
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(ggplot2)
library(gridExtra)

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine
Code
library(hrbrthemes)
library(psych)

Attaching package: 'psych'

The following objects are masked from 'package:ggplot2':

    %+%, alpha
Code
library(corrr)
library(lm.beta)
library (jtools)
library(broom.mixed)
library(car)
Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:psych':

    logit

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some
Code
# Read .csv file into a tibble
data <- read_csv("data/middle_school_music_mindsets_clean_partial.csv")
Rows: 503 Columns: 46
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (8): school, race, hispanic_latino, gender, currentSchoolMusic, school_...
dbl (38): survey_id, grade, schoolChoir, schoolInstrumental, schoolExtracurr...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Code
#reclassify character variables as factors

reclassFactor <- c("grade","school","race","hispanic_latino","gender","currentSchoolMusic","familyMusic")

data[reclassFactor] <- lapply(data[reclassFactor], as.factor)
Code
# Reorder the levels for currentSchoolMusic and set "No music" as the reference level
data$currentSchoolMusic <- relevel(data$currentSchoolMusic, ref = "No music")

Preliminary Analyses

Multicollinearity - VIF (Variance Inflation Factors)

I explored the predictors’ variance using Variance Inflation Factors (VIF) to identify potential collinearity. I followed the same steps for the three models (three different outcome variables).

Global MSC as outcome - VIF (Variance Inflation Factors)

Code
#VIF - All 17 predictors
car::vif(lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.524542  1        1.234724
active_engagement     1.491195  1        1.221145
musical_training      3.655350  1        1.911897
currentSchoolMusic    5.572385  3        1.331497
schoolChoir           1.740678  1        1.319348
schoolInstrumental    2.748262  1        1.657788
schoolExtracurricular 1.465815  1        1.210709
outsideChoir          1.237375  1        1.112374
outsideInstrumental   1.311194  1        1.145074
privateLesson         1.986178  1        1.409318
selfTaught            1.613293  1        1.270155
grade                 1.339242  2        1.075758
school                1.528026  2        1.111815
familyMusic           1.128720  1        1.062412
race                  1.448365  4        1.047393
hispanic_latino       1.188336  1        1.090108
gender                1.342277  3        1.050285

Based on these results, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.

Code
#VIF - After removing currentSchoolMusic
car::vif(lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.519985  1        1.232877
active_engagement     1.479603  1        1.216389
musical_training      3.500546  1        1.870975
schoolChoir           1.291890  1        1.136613
schoolInstrumental    1.684568  1        1.297909
schoolExtracurricular 1.441006  1        1.200419
outsideChoir          1.233339  1        1.110558
outsideInstrumental   1.302291  1        1.141180
privateLesson         1.898841  1        1.377984
selfTaught            1.607911  1        1.268034
grade                 1.225771  2        1.052210
school                1.173779  2        1.040870
familyMusic           1.111593  1        1.054321
race                  1.392029  4        1.042212
hispanic_latino       1.184485  1        1.088341
gender                1.296044  3        1.044167

After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.

Code
#VIF - After removing currentSchoolMusic and musical_training 
car::vif(lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.390065  1        1.179010
active_engagement     1.465463  1        1.210563
schoolChoir           1.269419  1        1.126685
schoolInstrumental    1.264796  1        1.124632
schoolExtracurricular 1.366288  1        1.168883
outsideChoir          1.233303  1        1.110542
outsideInstrumental   1.285356  1        1.133736
privateLesson         1.314733  1        1.146618
selfTaught            1.484374  1        1.218349
grade                 1.211958  2        1.049233
school                1.157202  2        1.037176
familyMusic           1.107391  1        1.052326
race                  1.367822  4        1.039929
hispanic_latino       1.184409  1        1.088306
gender                1.288264  3        1.043120

After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.

Ability MSC as coutcome - VIF (Variance Inflation Factors)

Code
#VIF - All 17 predictors
car::vif(lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.524542  1        1.234724
active_engagement     1.491195  1        1.221145
musical_training      3.655350  1        1.911897
currentSchoolMusic    5.572385  3        1.331497
schoolChoir           1.740678  1        1.319348
schoolInstrumental    2.748262  1        1.657788
schoolExtracurricular 1.465815  1        1.210709
outsideChoir          1.237375  1        1.112374
outsideInstrumental   1.311194  1        1.145074
privateLesson         1.986178  1        1.409318
selfTaught            1.613293  1        1.270155
grade                 1.339242  2        1.075758
school                1.528026  2        1.111815
familyMusic           1.128720  1        1.062412
race                  1.448365  4        1.047393
hispanic_latino       1.188336  1        1.090108
gender                1.342277  3        1.050285

Similar to the previous model, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.

Code
#VIF - After removing currentSchoolMusic
car::vif(lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.519985  1        1.232877
active_engagement     1.479603  1        1.216389
musical_training      3.500546  1        1.870975
schoolChoir           1.291890  1        1.136613
schoolInstrumental    1.684568  1        1.297909
schoolExtracurricular 1.441006  1        1.200419
outsideChoir          1.233339  1        1.110558
outsideInstrumental   1.302291  1        1.141180
privateLesson         1.898841  1        1.377984
selfTaught            1.607911  1        1.268034
grade                 1.225771  2        1.052210
school                1.173779  2        1.040870
familyMusic           1.111593  1        1.054321
race                  1.392029  4        1.042212
hispanic_latino       1.184485  1        1.088341
gender                1.296044  3        1.044167

After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.

Code
#VIF - After removing currentSchoolMusic and musical_training 
car::vif(lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.390065  1        1.179010
active_engagement     1.465463  1        1.210563
schoolChoir           1.269419  1        1.126685
schoolInstrumental    1.264796  1        1.124632
schoolExtracurricular 1.366288  1        1.168883
outsideChoir          1.233303  1        1.110542
outsideInstrumental   1.285356  1        1.133736
privateLesson         1.314733  1        1.146618
selfTaught            1.484374  1        1.218349
grade                 1.211958  2        1.049233
school                1.157202  2        1.037176
familyMusic           1.107391  1        1.052326
race                  1.367822  4        1.039929
hispanic_latino       1.184409  1        1.088306
gender                1.288264  3        1.043120

After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.

Non-Academic MSC as Outcome - VIF (Variance Inflation Factors)

Code
#VIF - All 17 predictors
car::vif(lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.524542  1        1.234724
active_engagement     1.491195  1        1.221145
musical_training      3.655350  1        1.911897
currentSchoolMusic    5.572385  3        1.331497
schoolChoir           1.740678  1        1.319348
schoolInstrumental    2.748262  1        1.657788
schoolExtracurricular 1.465815  1        1.210709
outsideChoir          1.237375  1        1.112374
outsideInstrumental   1.311194  1        1.145074
privateLesson         1.986178  1        1.409318
selfTaught            1.613293  1        1.270155
grade                 1.339242  2        1.075758
school                1.528026  2        1.111815
familyMusic           1.128720  1        1.062412
race                  1.448365  4        1.047393
hispanic_latino       1.188336  1        1.090108
gender                1.342277  3        1.050285

Consistent with the first two models, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.

Code
#VIF - After removing currentSchoolMusic
car::vif(lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.519985  1        1.232877
active_engagement     1.479603  1        1.216389
musical_training      3.500546  1        1.870975
schoolChoir           1.291890  1        1.136613
schoolInstrumental    1.684568  1        1.297909
schoolExtracurricular 1.441006  1        1.200419
outsideChoir          1.233339  1        1.110558
outsideInstrumental   1.302291  1        1.141180
privateLesson         1.898841  1        1.377984
selfTaught            1.607911  1        1.268034
grade                 1.225771  2        1.052210
school                1.173779  2        1.040870
familyMusic           1.111593  1        1.054321
race                  1.392029  4        1.042212
hispanic_latino       1.184485  1        1.088341
gender                1.296044  3        1.044167

After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.

Code
#VIF - After removing currentSchoolMusic and musical_training 
car::vif(lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data))
                          GVIF Df GVIF^(1/(2*Df))
growth_mindset        1.390065  1        1.179010
active_engagement     1.465463  1        1.210563
schoolChoir           1.269419  1        1.126685
schoolInstrumental    1.264796  1        1.124632
schoolExtracurricular 1.366288  1        1.168883
outsideChoir          1.233303  1        1.110542
outsideInstrumental   1.285356  1        1.133736
privateLesson         1.314733  1        1.146618
selfTaught            1.484374  1        1.218349
grade                 1.211958  2        1.049233
school                1.157202  2        1.037176
familyMusic           1.107391  1        1.052326
race                  1.367822  4        1.039929
hispanic_latino       1.184409  1        1.088306
gender                1.288264  3        1.043120

After removing currentSchoolMusic and musical_training, there seem to be no collinearity between the predictors. Hence,I will not include those variables in the model.

Assumptions

Before building the linear regression models, I examined assumptions of linearity, homoscedasticity, and other assumptions. I also examine if using the music self-concept composite score as an outcome would yield different results compared to the musical ability/music self-concept subscale.

Code
#Global MSC Final Model
final_model_msc_full <- lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)

#Ability MSC Final Model
final_model_msc_ability <- lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)

#No Ability MSC Final Model
final_model_msc_no_ability <- lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)

Normality

Q-Q Plot

Code
qqPlot(final_model_msc_full)

[1]  47 336
Code
qqPlot(final_model_msc_ability)

[1] 104 307
Code
qqPlot(final_model_msc_no_ability)

[1] 112 211

Kolmogorov-Smirnov

Code
# Global MSC
ks.test(residuals(lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  residuals(lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
D = 0.39939, p-value < 2.2e-16
alternative hypothesis: two-sided
Code
#Ability MSC

ks.test(residuals(lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  residuals(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
D = 0.1943, p-value < 2.2e-16
alternative hypothesis: two-sided
Code
#No Ability MSC
ks.test(residuals(lm(formula =  msc_no_ability ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  residuals(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
D = 0.39233, p-value < 2.2e-16
alternative hypothesis: two-sided

Shapiro-Wilk

Code
# Global MSC
shapiro.test(residuals(lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)))

    Shapiro-Wilk normality test

data:  residuals(lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
W = 0.9932, p-value = 0.02254
Code
#Ability MSC

shapiro.test(residuals(lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)))

    Shapiro-Wilk normality test

data:  residuals(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
W = 0.99464, p-value = 0.07616
Code
#No Ability MSC
shapiro.test(residuals(lm(formula =  msc_no_ability ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)))

    Shapiro-Wilk normality test

data:  residuals(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))
W = 0.99145, p-value = 0.005353

Kolmogorov-Smirnov and Shapiro- Wilk tests showed that the data follows a nonnormal distribution. Individual examination of each variables’ skewnewss and kurtosis identify three variables, but the distribution was not extreme.

Stepwise Multiple Linear Regression

I conducted a series of linear regression models to examine the relationships between multiple predictors and each of the three outcome variables, global MSC, ability MSC, and non-academic MSC. The aim was to address the research questions gradually, facilitating the interpretability of the findings by assessing the influence of successive categories of predictors on the models’ improvement. I conducted the same steps and used the same predictors for each outcome variable.

On Step 1, I introduced mindset of music ability as the sole predictor.

To investigate the impact of various types of music engagement on music self-concept, I introduced two sets of variables in the following steps. Step 2 included active engagement, whereas Step 3 comprised variables related to duration of engagement in music activities in and out of school.

In the following and final step (Step 4), I introduced background variables like grade, school, race, ethnicity, gender, and family music background.

GLOBAL MSC - Stepwise Multiple Linear Regression

Intercept

Code
#Intercept
lm_msc_full_step_0 <- lm(formula = msc_full ~ 1, data = data)

summary(lm_msc_full_step_0)

Call:
lm(formula = msc_full ~ 1, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.748  -8.748   0.252   9.252  32.252 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  75.7475     0.6215   121.9   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.94 on 502 degrees of freedom

Step 1

Code
#Step 1
lm_msc_full_step_1 <- lm(formula = msc_full ~ growth_mindset,  data = data)

summary_lm_msc_full_step_1 <- summary(lm_msc_full_step_1)

summary_lm_msc_full_step_1

Call:
lm(formula = msc_full ~ growth_mindset, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-41.934  -7.934  -0.037   8.066  31.859 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     42.5431     2.5017   17.00   <2e-16 ***
growth_mindset   2.6494     0.1951   13.58   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.93 on 501 degrees of freedom
Multiple R-squared:  0.2691,    Adjusted R-squared:  0.2677 
F-statistic: 184.5 on 1 and 501 DF,  p-value: < 2.2e-16

Mindset of music ability significantly predicted global MSC, explaining 27% of its variance (adjusted R2 = 0.2677, p < .01).

Step 2

Code
#Step 2
lm_msc_full_step_2 <- lm(formula = msc_full ~ growth_mindset  +
                   active_engagement,  data = data)

summary_lm_msc_full_step_2 <- summary(lm_msc_full_step_2)

summary_lm_msc_full_step_2

Call:
lm(formula = msc_full ~ growth_mindset + active_engagement, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-36.328  -5.335   0.328   5.410  26.347 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       17.31352    2.20194   7.863 2.33e-14 ***
growth_mindset     1.79947    0.14887  12.087  < 2e-16 ***
active_engagement  1.56372    0.07537  20.748  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.752 on 500 degrees of freedom
Multiple R-squared:  0.6073,    Adjusted R-squared:  0.6057 
F-statistic: 386.6 on 2 and 500 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_full_1_2 <- anova(lm_msc_full_step_1, lm_msc_full_step_2)

f_test_result_lm_msc_full_1_2
Analysis of Variance Table

Model 1: msc_full ~ growth_mindset
Model 2: msc_full ~ growth_mindset + active_engagement
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    501 71272                                  
2    500 38299  1     32974 430.48 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In step 2, active engagement emerged as a significant predictor, collectively contributing to 62% of the variance in MSC (adjusted R2 = 0.6057, p < .001 ). The addition of these variables significantly improved the model, F = 430.48, p < .001).

Step 3

Code
#Step 3
lm_msc_full_step_3 <- lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + 
                   outsideChoir + outsideInstrumental + 
                   privateLesson + selfTaught,  data = data)

summary_lm_msc_full_step_3 <- summary(lm_msc_full_step_3)

summary_lm_msc_full_step_3

Call:
lm(formula = msc_full ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.015  -5.394   0.183   5.203  27.355 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           21.05592    2.20238   9.561  < 2e-16 ***
growth_mindset         1.57414    0.15446  10.191  < 2e-16 ***
active_engagement      1.36259    0.07879  17.294  < 2e-16 ***
schoolChoir            0.10021    0.22212   0.451   0.6521    
schoolInstrumental    -0.36064    0.24906  -1.448   0.1482    
schoolExtracurricular  0.60838    0.28759   2.115   0.0349 *  
outsideChoir           0.57266    0.42542   1.346   0.1789    
outsideInstrumental   -0.36545    0.48373  -0.755   0.4503    
privateLesson         -0.06748    0.20658  -0.327   0.7441    
selfTaught             1.26614    0.22971   5.512 5.73e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.358 on 493 degrees of freedom
Multiple R-squared:  0.6468,    Adjusted R-squared:  0.6404 
F-statistic: 100.3 on 9 and 493 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_full_2_3 <- anova(lm_msc_full_step_2, lm_msc_full_step_3)

f_test_result_lm_msc_full_2_3
Analysis of Variance Table

Model 1: msc_full ~ growth_mindset + active_engagement
Model 2: msc_full ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    500 38299                                  
2    493 34441  7    3857.8 7.8888 4.437e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The refined model explained 65% of the variance in music self-concept (adjusted R2 = 0.6404, p < .01), representing a significant fit improvement, F = 7.89, p < .001).

Step 4

Code
#Step 4
lm_msc_full_step_4 <- lm(formula = msc_full ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)

summary_lm_msc_full_step_4 <- summary(lm.beta(lm_msc_full_step_4))

summary_lm_msc_full_step_4

Call:
lm(formula = msc_full ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught + 
    grade + school + familyMusic + race + hispanic_latino + gender, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.838  -5.277   0.388   4.854  25.683 

Coefficients:
                               Estimate Standardized Std. Error t value
(Intercept)                   24.364574           NA   3.198147   7.618
growth_mindset                 1.579534     0.309296   0.160000   9.872
active_engagement              1.303978     0.504367   0.083169  15.679
schoolChoir                    0.151462     0.019663   0.230619   0.657
schoolInstrumental            -0.202408    -0.023435   0.258114  -0.784
schoolExtracurricular          0.576783     0.061976   0.289076   1.995
outsideChoir                   0.443235     0.030298   0.431720   1.027
outsideInstrumental           -0.313472    -0.019189   0.492157  -0.637
privateLesson                 -0.065704    -0.009427   0.212359  -0.309
selfTaught                     1.255874     0.175800   0.231285   5.430
grade7                        -1.160866    -0.038698   0.967888  -1.199
grade8                        -1.421237    -0.049382   0.966866  -1.470
schoolB                       -0.718357    -0.025781   0.950530  -0.756
schoolC                        0.397597     0.012653   1.074343   0.370
familyMusicYes                -0.082238    -0.002943   0.781334  -0.105
raceBlack or African American  2.524285     0.050152   2.426954   1.040
raceSome Other Race            1.090157     0.015299   2.769650   0.394
raceTwo or More Races          1.242023     0.022994   2.515745   0.494
raceWhite                     -0.208615    -0.006311   2.062744  -0.101
hispanic_latinoYes            -0.517686    -0.007925   1.889096  -0.274
genderMale                    -2.475762    -0.088734   0.817231  -3.029
genderNon-Binary              -2.947811    -0.018804   4.280852  -0.689
genderPrefer not to say       -5.901356    -0.053024   3.029488  -1.948
                              Pr(>|t|)    
(Intercept)                   1.38e-13 ***
growth_mindset                 < 2e-16 ***
active_engagement              < 2e-16 ***
schoolChoir                    0.51165    
schoolInstrumental             0.43332    
schoolExtracurricular          0.04658 *  
outsideChoir                   0.30509    
outsideInstrumental            0.52447    
privateLesson                  0.75715    
selfTaught                    8.97e-08 ***
grade7                         0.23097    
grade8                         0.14223    
schoolB                        0.45017    
schoolC                        0.71148    
familyMusicYes                 0.91622    
raceBlack or African American  0.29882    
raceSome Other Race            0.69405    
raceTwo or More Races          0.62174    
raceWhite                      0.91949    
hispanic_latinoYes             0.78417    
genderMale                     0.00258 ** 
genderNon-Binary               0.49141    
genderPrefer not to say        0.05200 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.298 on 480 degrees of freedom
Multiple R-squared:  0.661, Adjusted R-squared:  0.6455 
F-statistic: 42.55 on 22 and 480 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_full_3_4 <- anova(lm_msc_full_step_3, lm_msc_full_step_4)

f_test_result_lm_msc_full_3_4
Analysis of Variance Table

Model 1: msc_full ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
Model 2: msc_full ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1    493 34441                              
2    480 33054 13    1386.7 1.5491 0.09643 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In the new model, gender was a significant predictor of global music self-concept, with male students differing significantly from female students (reference group). This extended model explained 68% of the variation in music self-concept (adjusted R2 = 0.6455, p = .049), showing a significantly superior fit compared to the previous model (F = 1.55, p < .01).

Results Table

Code
#RESULTS TABLE Global MSC

# Extract coefficients and statistics
results_lm_msc_full <- bind_rows(broom::tidy(lm.beta(lm_msc_full_step_4), conf.int = TRUE)
)

# Round output to two decimal points
results_lm_msc_full$estimate <- round(results_lm_msc_full$estimate,2)

results_lm_msc_full$std_estimate <- round(results_lm_msc_full$std_estimate,3)

results_lm_msc_full$std.error <- round(results_lm_msc_full$std.error,2)

results_lm_msc_full$statistic <- round(results_lm_msc_full$statistic,2)

results_lm_msc_full$p.value <- round(results_lm_msc_full$p.value,3)

results_lm_msc_full$conf.low <- round(results_lm_msc_full$conf.low,2)

results_lm_msc_full$conf.high <- round(results_lm_msc_full$conf.high,2)


# Combine conf.low and conf.high into a single column called "95% CI" and change column names
results_lm_msc_full <- results_lm_msc_full %>%
  mutate("95% CI" = paste0("[", conf.low, ", ", conf.high, "]")) %>%
  select(-conf.low, -conf.high) %>%   # Remove the original conf.low and conf.high columns
  rename(b = estimate, 
         β = std_estimate,
           SE = std.error,
         t = statistic,
         p = p.value)

# Change order of columns
results_lm_msc_full <- results_lm_msc_full %>% 
  select(term, b, SE, β, t, "95% CI", p)


# Change output for p values undeer .001
results_lm_msc_full <- results_lm_msc_full %>% 
  mutate(p = ifelse(p < 0.001, "<.001", p))

# Change predictors' names
results_lm_msc_full <- results_lm_msc_full %>% 
  mutate(term = case_match(term, 
                      "(Intercept)" ~ "(Intercept)",
                      "growth_mindset" ~ "Growth mindset", 
                      "active_engagement" ~ "Active Engagement",
                      "schoolChoir" ~ "School choir",
                      "schoolInstrumental" ~ "School instrumental",
                      "schoolExtracurricular" ~ "School music club",
                      "outsideChoir" ~ "Choir (outside school)",
                      "outsideInstrumental" ~ "instrumental (outside school)",
                      "privateLesson" ~ "Private lesson",
                      "selfTaught" ~ "Self-teaching",
                      "grade7" ~ "Grade 7",
                      "grade8" ~ "Grade 8",
                      "schoolB" ~ "School B",
                      "schoolC" ~ "School C",
                      "familyMusicYes" ~ "Family music engagement",
                      "raceBlack or African American" ~ "Black or African American",
                      "raceSome Other Race" ~ "Other race",
                      "raceTwo or More Races" ~ "Two or more races",
                      "raceWhite" ~ "White",
                      "hispanic_latinoYes" ~ "Ethnicity",
                      "genderMale" ~ "Male",
                      "genderNon-Binary" ~ "Non-Binary",
                      "genderPrefer not to say" ~ "Gender not disclosed")
                      )

print(results_lm_msc_full, n= 100)
# A tibble: 23 × 7
   term                              b    SE      β     t `95% CI`      p    
   <chr>                         <dbl> <dbl>  <dbl> <dbl> <chr>         <chr>
 1 (Intercept)                   24.4   3.2  NA      7.62 [NA, NA]      <.001
 2 Growth mindset                 1.58  0.16  0.309  9.87 [-0.01, 0.62] <.001
 3 Active Engagement              1.3   0.08  0.504 15.7  [0.34, 0.67]  <.001
 4 School choir                   0.15  0.23  0.02   0.66 [-0.43, 0.47] 0.512
 5 School instrumental           -0.2   0.26 -0.023 -0.78 [-0.53, 0.48] 0.433
 6 School music club              0.58  0.29  0.062  2    [-0.51, 0.63] 0.047
 7 Choir (outside school)         0.44  0.43  0.03   1.03 [-0.82, 0.88] 0.305
 8 instrumental (outside school) -0.31  0.49 -0.019 -0.64 [-0.99, 0.95] 0.524
 9 Private lesson                -0.07  0.21 -0.009 -0.31 [-0.43, 0.41] 0.757
10 Self-teaching                  1.26  0.23  0.176  5.43 [-0.28, 0.63] <.001
11 Grade 7                       -1.16  0.97 -0.039 -1.2  [-1.94, 1.86] 0.231
12 Grade 8                       -1.42  0.97 -0.049 -1.47 [-1.95, 1.85] 0.142
13 School B                      -0.72  0.95 -0.026 -0.76 [-1.89, 1.84] 0.45 
14 School C                       0.4   1.07  0.013  0.37 [-2.1, 2.12]  0.711
15 Family music engagement       -0.08  0.78 -0.003 -0.11 [-1.54, 1.53] 0.916
16 Black or African American      2.52  2.43  0.05   1.04 [-4.72, 4.82] 0.299
17 Other race                     1.09  2.77  0.015  0.39 [-5.43, 5.46] 0.694
18 Two or more races              1.24  2.52  0.023  0.49 [-4.92, 4.97] 0.622
19 White                         -0.21  2.06 -0.006 -0.1  [-4.06, 4.05] 0.919
20 Ethnicity                     -0.52  1.89 -0.008 -0.27 [-3.72, 3.7]  0.784
21 Male                          -2.48  0.82 -0.089 -3.03 [-1.69, 1.52] 0.003
22 Non-Binary                    -2.95  4.28 -0.019 -0.69 [-8.43, 8.39] 0.491
23 Gender not disclosed          -5.9   3.03 -0.053 -1.95 [-6.01, 5.9]  0.052

Overall Significance

Code
#Calculate overall significance/fit

f_test_result_lm_msc_full_0_4 <- anova(lm_msc_full_step_0, lm_msc_full_step_4)

f_test_result_lm_msc_full_0_4
Analysis of Variance Table

Model 1: msc_full ~ 1
Model 2: msc_full ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    502 97519                                  
2    480 33054 22     64465 42.551 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The overall significance of the final model was confirmed through the F-test (F = 42.55, p < .01).

Post Hoc Analyses

Equality of coefficients (Wald Test)

Although several predictors showed significant effect on the outcome, their coefficient may or may not differ from each other significantly. To examine if the differences between the coefficients (as observed in values of b and β ) are significant, I visually inspected those coefficients.

Code
#Visually compare coefficients
lm_full <- lm(formula = msc_full ~ growth_mindset + active_engagement + selfTaught + gender, data = data)

jtools::plot_summs(lm_full)


Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.

Code
wt_mindset_active <- linearHypothesis(lm_full, "growth_mindset  = active_engagement ")

wt_mindset_active

Linear hypothesis test:
growth_mindset - active_engagement = 0

Model 1: restricted model
Model 2: msc_full ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1    497 34487                              
2    496 34276  1    210.59 3.0474 0.08148 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_mindset_selfTaught <- linearHypothesis(lm_full, "growth_mindset  = selfTaught ")

wt_mindset_selfTaught

Linear hypothesis test:
growth_mindset - selfTaught = 0

Model 1: restricted model
Model 2: msc_full ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    497 34312                           
2    496 34276  1    35.847 0.5187 0.4717

Code
wt_active_selfTaught <- linearHypothesis(lm_full, "active_engagement  = selfTaught ")

wt_active_selfTaught

Linear hypothesis test:
active_engagement - selfTaught = 0

Model 1: restricted model
Model 2: msc_full ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    497 34290                           
2    496 34276  1     13.72 0.1985 0.6561

Wald tests results indicated that the difference in coefficients between mindset of music ability and active engagement, mindset of music abilty and self-teaching, and active engagement and self-teaching were not statistically significant, suggesting those predictors had the same effect on students’ global music self-concept.

Joint Significance (Wald Test)

Code
wt_gender <- lmtest::waldtest(lm_full, term = "gender")

print(wt_gender)
Wald test

Model 1: msc_full ~ growth_mindset + active_engagement + selfTaught + 
    gender
Model 2: msc_full ~ growth_mindset + active_engagement + selfTaught
  Res.Df Df      F  Pr(>F)   
1    496                     
2    499 -3 4.2503 0.00558 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I employed the Wald test once more to delve deeper into the joint significance of the gender variable. The results showed a significant difference between the model incorporating gender and the one excluding it, F = 4.2502945, p = < .01. This finding indicates that gender stands as a significant predictor of music self-concept, underscoring its importance in the predictive model.

Scatter Plots

Code
#scatter plots with fitted regression lines
plot1_full <- ggplot(data, aes(x = growth_mindset_z,  y = msc_full_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Mindset of Music Ability",
       y = "Overall") +
  theme_ipsum()

plot2_full <- ggplot(data, aes(x = active_engagement_z,  y = msc_full_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Active Engagement",
       y = "Global Music Self-Concept") +
  theme_ipsum()


plot3_full <- ggplot(data, aes(x = selfTaught_z,  y = msc_full_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Self-Teaching",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot4_full <- ggplot(data, aes(x = schoolExtracurricular_z,  y = msc_full_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Extracurricular School Music Clubs",
       y = "Global Music Self-Concept") +
  theme_ipsum()

#Combine plots into a single display (grid layout)
grid.arrange(plot1_full, plot2_full, plot3_full, plot4_full, ncol = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'


ACADEMIC MSC (Music Ability Only) - Stepwise Multiple Linear Regression

Intercept

Code
#Intercept
lm_msc_ability_step_0 <- lm(formula = musical_ability_msc ~ 1, data = data)

summary(lm_msc_ability_step_0)

Call:
lm(formula = musical_ability_msc ~ 1, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-7.173 -3.173 -0.173  2.827  7.827 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.1730     0.1603   75.92   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.596 on 502 degrees of freedom

Step 1

Code
#Step 1
lm_msc_ability_step_1 <- lm(formula = musical_ability_msc ~ growth_mindset,  data = data)

summary_lm_msc_ability_step_1 <- summary(lm_msc_ability_step_1)

summary_lm_msc_ability_step_1

Call:
lm(formula = musical_ability_msc ~ growth_mindset, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8255 -1.8255  0.1745  2.1745  9.3550 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     2.58486    0.61474   4.205  3.1e-05 ***
growth_mindset  0.76504    0.04793  15.962  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.931 on 501 degrees of freedom
Multiple R-squared:  0.3371,    Adjusted R-squared:  0.3358 
F-statistic: 254.8 on 1 and 501 DF,  p-value: < 2.2e-16

Step 2

Code
#Step 2
lm_msc_ability_step_2 <- lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement,  data = data)

summary_lm_msc_ability_step_2 <- summary(lm_msc_ability_step_2)

summary_lm_msc_ability_step_2

Call:
lm(formula = musical_ability_msc ~ growth_mindset + active_engagement, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.5612 -1.8720  0.0405  1.8709  9.4115 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.37141    0.69821  -0.532    0.595    
growth_mindset     0.66545    0.04721  14.097  < 2e-16 ***
active_engagement  0.18323    0.02390   7.667 9.24e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.775 on 500 degrees of freedom
Multiple R-squared:  0.4068,    Adjusted R-squared:  0.4045 
F-statistic: 171.5 on 2 and 500 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_ability_1_2 <- anova(lm_msc_ability_step_1, lm_msc_ability_step_2)

f_test_result_lm_msc_ability_1_2
Analysis of Variance Table

Model 1: musical_ability_msc ~ growth_mindset
Model 2: musical_ability_msc ~ growth_mindset + active_engagement
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    501 4303.4                                  
2    500 3850.7  1    452.72 58.784 9.239e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 3

Code
#Step 3
lm_msc_ability_step_3 <- lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + 
                   outsideChoir + outsideInstrumental + 
                   privateLesson + selfTaught,  data = data)

summary_lm_msc_ability_step_3 <- summary(lm_msc_ability_step_3)

summary_lm_msc_ability_step_3

Call:
lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.4186 -1.6874  0.0977  1.6702  8.4095 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            1.28572    0.64846   1.983 0.047953 *  
growth_mindset         0.48158    0.04548  10.589  < 2e-16 ***
active_engagement      0.09627    0.02320   4.150 3.92e-05 ***
schoolChoir            0.01977    0.06540   0.302 0.762595    
schoolInstrumental     0.27304    0.07333   3.723 0.000219 ***
schoolExtracurricular  0.42300    0.08468   4.996 8.16e-07 ***
outsideChoir           0.22050    0.12526   1.760 0.078962 .  
outsideInstrumental   -0.34493    0.14243  -2.422 0.015804 *  
privateLesson          0.16628    0.06082   2.734 0.006486 ** 
selfTaught             0.44399    0.06763   6.565 1.33e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.461 on 493 degrees of freedom
Multiple R-squared:  0.5401,    Adjusted R-squared:  0.5317 
F-statistic: 64.33 on 9 and 493 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_ability_2_3 <- anova(lm_msc_ability_step_2, lm_msc_ability_step_3)

f_test_result_lm_msc_ability_2_3
Analysis of Variance Table

Model 1: musical_ability_msc ~ growth_mindset + active_engagement
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    500 3850.7                                  
2    493 2985.8  7    864.95 20.402 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 4

Code
#Step 4
lm_msc_ability_step_4 <- lm(formula = musical_ability_msc ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  +
                   grade + school + familyMusic + race + hispanic_latino + gender, data = data)

summary_lm_msc_ability_step_4 <- summary(lm.beta(lm_msc_ability_step_4))

summary_lm_msc_ability_step_4

Call:
lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught + 
    grade + school + familyMusic + race + hispanic_latino + gender, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5097 -1.3775  0.1049  1.5720  9.1139 

Coefficients:
                               Estimate Standardized Std. Error t value
(Intercept)                    1.977131           NA   0.907403   2.179
growth_mindset                 0.407622     0.309357   0.045396   8.979
active_engagement              0.127718     0.191463   0.023597   5.412
schoolChoir                    0.088781     0.044672   0.065433   1.357
schoolInstrumental             0.295989     0.132825   0.073234   4.042
schoolExtracurricular          0.371470     0.154699   0.082019   4.529
outsideChoir                   0.187682     0.049723   0.122491   1.532
outsideInstrumental           -0.414830    -0.098420   0.139639  -2.971
privateLesson                  0.192790     0.107211   0.060252   3.200
selfTaught                     0.440988     0.239252   0.065622   6.720
grade7                        -0.547205    -0.070700   0.274616  -1.993
grade8                        -0.945464    -0.127323   0.274327  -3.446
schoolB                       -0.278372    -0.038721   0.269691  -1.032
schoolC                        0.850315     0.104881   0.304821   2.790
familyMusicYes                 0.277356     0.038473   0.221686   1.251
raceBlack or African American  0.178241     0.013725   0.688594   0.259
raceSome Other Race           -0.250633    -0.013632   0.785826  -0.319
raceTwo or More Races         -0.022927    -0.001645   0.713786  -0.032
raceWhite                     -0.677802    -0.079472   0.585258  -1.158
hispanic_latinoYes            -0.908193    -0.053886   0.535989  -1.694
genderMale                     0.711032     0.098771   0.231871   3.066
genderNon-Binary               0.551157     0.013626   1.214596   0.454
genderPrefer not to say       -1.449794    -0.050487   0.859549  -1.687
                              Pr(>|t|)    
(Intercept)                   0.029825 *  
growth_mindset                 < 2e-16 ***
active_engagement             9.84e-08 ***
schoolChoir                   0.175477    
schoolInstrumental            6.18e-05 ***
schoolExtracurricular         7.48e-06 ***
outsideChoir                  0.126129    
outsideInstrumental           0.003120 ** 
privateLesson                 0.001467 ** 
selfTaught                    5.16e-11 ***
grade7                        0.046869 *  
grade8                        0.000618 ***
schoolB                       0.302504    
schoolC                       0.005488 ** 
familyMusicYes                0.211500    
raceBlack or African American 0.795864    
raceSome Other Race           0.749909    
raceTwo or More Races         0.974389    
raceWhite                     0.247389    
hispanic_latinoYes            0.090833 .  
genderMale                    0.002288 ** 
genderNon-Binary              0.650193    
genderPrefer not to say       0.092312 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.354 on 480 degrees of freedom
Multiple R-squared:  0.5901,    Adjusted R-squared:  0.5713 
F-statistic: 31.41 on 22 and 480 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_ability_3_4 <- anova(lm_msc_ability_step_3, lm_msc_ability_step_4)

f_test_result_lm_msc_ability_3_4
Analysis of Variance Table

Model 1: musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    493 2985.8                                  
2    480 2660.9 13    324.86 4.5079 2.787e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results Table

Code
#RESULTS TABLE: Ability MSC

# Extract coefficients and statistics
results_lm_msc_ability <- bind_rows(broom::tidy(lm.beta(lm_msc_ability_step_4), conf.int = TRUE)
)


# Round output to two decimal points
results_lm_msc_ability$estimate <- round(results_lm_msc_ability$estimate,2)

results_lm_msc_ability$std_estimate <- round(results_lm_msc_ability$std_estimate,2)

results_lm_msc_ability$std.error <- round(results_lm_msc_ability$std.error,2)

results_lm_msc_ability$statistic <- round(results_lm_msc_ability$statistic,2)

results_lm_msc_ability$p.value <- round(results_lm_msc_ability$p.value,3)

results_lm_msc_ability$conf.low <- round(results_lm_msc_ability$conf.low,2)

results_lm_msc_ability$conf.high <- round(results_lm_msc_ability$conf.high,2)


# Combine conf.low and conf.high into a single column called "95% CI" and change column names
results_lm_msc_ability <- results_lm_msc_ability %>%
  mutate("95% CI" = paste0("[", conf.low, ", ", conf.high, "]")) %>%
  select(-conf.low, -conf.high) %>%   # Remove the original conf.low and conf.high columns
  rename(b = estimate, 
         β = std_estimate,
           SE = std.error,
         t = statistic,
         p = p.value)

# Change order of columns
results_lm_msc_ability <- results_lm_msc_ability %>% 
  select(term, b, SE, β, t, "95% CI", p)


# Change output for p values undeer .001
results_lm_msc_ability <- results_lm_msc_ability %>% 
  mutate(p = ifelse(p < 0.001, "<.001", p))

# Change predictors' names
results_lm_msc_ability <- results_lm_msc_ability %>% 
  mutate(term = 
           case_match(term, 
                      "(Intercept)" ~ "(Intercept)",
                      "growth_mindset" ~ "Growth mindset", 
                      "active_engagement" ~ "Active Engagement",
                      "schoolChoir" ~ "School choir",
                      "schoolInstrumental" ~ "School instrumental",
                      "schoolExtracurricular" ~ "School music club",
                      "outsideChoir" ~ "Choir (outside school)",
                      "outsideInstrumental" ~ "instrumental (outside school)",
                      "privateLesson" ~ "Private lesson",
                      "selfTaught" ~ "Self-teaching",
                      "grade7" ~ "Grade 7",
                      "grade8" ~ "Grade 8",
                      "schoolB" ~ "School B",
                      "schoolC" ~ "School C",
                      "familyMusicYes" ~ "Family music engagement",
                      "raceBlack or African American" ~ "Black or African American",
                      "raceSome Other Race" ~ "Other race",
                      "raceTwo or More Races" ~ "Two or more races",
                      "raceWhite" ~ "White",
                      "hispanic_latinoYes" ~ "Ethnicity",
                      "genderMale" ~ "Male",
                      "genderNon-Binary" ~ "Non-Binary",
                      "genderPrefer not to say" ~ "Gender not disclosed")
                      )

print(results_lm_msc_ability, n= 100)
# A tibble: 23 × 7
   term                              b    SE     β     t `95% CI`      p    
   <chr>                         <dbl> <dbl> <dbl> <dbl> <chr>         <chr>
 1 (Intercept)                    1.98  0.91 NA     2.18 [NA, NA]      0.03 
 2 Growth mindset                 0.41  0.05  0.31  8.98 [0.22, 0.4]   <.001
 3 Active Engagement              0.13  0.02  0.19  5.41 [0.15, 0.24]  <.001
 4 School choir                   0.09  0.07  0.04  1.36 [-0.08, 0.17] 0.175
 5 School instrumental            0.3   0.07  0.13  4.04 [-0.01, 0.28] <.001
 6 School music club              0.37  0.08  0.15  4.53 [-0.01, 0.32] <.001
 7 Choir (outside school)         0.19  0.12  0.05  1.53 [-0.19, 0.29] 0.126
 8 instrumental (outside school) -0.41  0.14 -0.1  -2.97 [-0.37, 0.18] 0.003
 9 Private lesson                 0.19  0.06  0.11  3.2  [-0.01, 0.23] 0.001
10 Self-teaching                  0.44  0.07  0.24  6.72 [0.11, 0.37]  <.001
11 Grade 7                       -0.55  0.27 -0.07 -1.99 [-0.61, 0.47] 0.047
12 Grade 8                       -0.95  0.27 -0.13 -3.45 [-0.67, 0.41] 0.001
13 School B                      -0.28  0.27 -0.04 -1.03 [-0.57, 0.49] 0.303
14 School C                       0.85  0.3   0.1   2.79 [-0.49, 0.7]  0.005
15 Family music engagement        0.28  0.22  0.04  1.25 [-0.4, 0.47]  0.211
16 Black or African American      0.18  0.69  0.01  0.26 [-1.34, 1.37] 0.796
17 Other race                    -0.25  0.79 -0.01 -0.32 [-1.56, 1.53] 0.75 
18 Two or more races             -0.02  0.71  0    -0.03 [-1.4, 1.4]   0.974
19 White                         -0.68  0.59 -0.08 -1.16 [-1.23, 1.07] 0.247
20 Ethnicity                     -0.91  0.54 -0.05 -1.69 [-1.11, 1]    0.091
21 Male                           0.71  0.23  0.1   3.07 [-0.36, 0.55] 0.002
22 Non-Binary                     0.55  1.21  0.01  0.45 [-2.37, 2.4]  0.65 
23 Gender not disclosed          -1.45  0.86 -0.05 -1.69 [-1.74, 1.64] 0.092

Overall Significance

Code
#overall significance/fit
f_test_result_lm_msc_ability_0_4 <- anova(lm_msc_ability_step_0, lm_msc_ability_step_4)

f_test_result_lm_msc_ability_0_4
Analysis of Variance Table

Model 1: musical_ability_msc ~ 1
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    502 6492.0                                  
2    480 2660.9 22      3831 31.413 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post Hoc Analyses

Equality of coefficients (Wald Test)

Code
#Visually compare coefficients
library (jtools)
library(broom.mixed)

lm_ability <- lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + schoolExtracurricular + privateLesson + selfTaught + grade + school + gender, data = data)

jtools::plot_summs(lm_ability)


Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, school instrumental ensemble, extracurricular school music club, private lesson, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.

Code
library(car)

wt_ability_mindset_active <- linearHypothesis(lm_ability, "growth_mindset  = active_engagement ")

wt_ability_mindset_active

Linear hypothesis test:
growth_mindset - active_engagement = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    490 3008.9                                  
2    489 2815.1  1    193.83 33.669 1.175e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_mindset_schoolInstrumental <- linearHypothesis(lm_ability, "growth_mindset  = schoolInstrumental ")

wt_ability_mindset_schoolInstrumental

Linear hypothesis test:
growth_mindset - schoolInstrumental = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    490 2844.2                              
2    489 2815.1  1    29.088 5.0527 0.02503 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_mindset_schoolExtracurricular <- linearHypothesis(lm_ability, "growth_mindset  = schoolExtracurricular ")

wt_ability_mindset_schoolExtracurricular

Linear hypothesis test:
growth_mindset - schoolExtracurricular = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2819.5                           
2    489 2815.1  1     4.406 0.7653 0.3821

Code
wt_ability_mindset_privateLesson <- linearHypothesis(lm_ability, "growth_mindset  = privateLesson ")

wt_ability_mindset_privateLesson

Linear hypothesis test:
growth_mindset - privateLesson = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    490 2871.5                                
2    489 2815.1  1    56.428 9.8018 0.001848 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_mindset_selfTaught <- linearHypothesis(lm_ability, "growth_mindset  = selfTaught ")

wt_ability_mindset_selfTaught

Linear hypothesis test:
growth_mindset - selfTaught = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2815.5                           
2    489 2815.1  1   0.39929 0.0694 0.7924

Code
wt_ability_active_schoolInstrumental <- linearHypothesis(lm_ability, "active_engagement  = schoolInstrumental ")

wt_ability_active_schoolInstrumental

Linear hypothesis test:
active_engagement - schoolInstrumental = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1    490 2828.3                          
2    489 2815.1  1    13.218 2.296 0.1304

Code
wt_ability_active_schoolExtracurricular <- linearHypothesis(lm_ability, "active_engagement  = schoolExtracurricular ")

wt_ability_active_schoolExtracurricular

Linear hypothesis test:
active_engagement - schoolExtracurricular = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    490 2857.5                                
2    489 2815.1  1     42.42 7.3686 0.006872 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_active_privateLesson <- linearHypothesis(lm_ability, "active_engagement  = privateLesson ")

wt_ability_active_privateLesson

Linear hypothesis test:
active_engagement - privateLesson = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2820.6                           
2    489 2815.1  1    5.4757 0.9512 0.3299

Code
wt_ability_active_selfTaught <- linearHypothesis(lm_ability, "active_engagement  = selfTaught ")

wt_ability_active_selfTaught

Linear hypothesis test:
active_engagement - selfTaught = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    490 2935.7                                  
2    489 2815.1  1    120.62 20.952 5.981e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_schoolInstrumental_schoolExtracurricular <- linearHypothesis(lm_ability, "schoolInstrumental  = schoolExtracurricular ")

wt_ability_schoolInstrumental_schoolExtracurricular

Linear hypothesis test:
schoolInstrumental - schoolExtracurricular = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2822.0                           
2    489 2815.1  1    6.8768 1.1945  0.275

Code
wt_ability_schoolInstrumental_privateLesson <- linearHypothesis(lm_ability, "schoolInstrumental  = privateLesson ")

wt_ability_schoolInstrumental_privateLesson

Linear hypothesis test:
schoolInstrumental - privateLesson = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2816.7                           
2    489 2815.1  1    1.5696 0.2726 0.6018

Code
wt_ability_schoolInstrumental_selfTaught <- linearHypothesis(lm_ability, "schoolInstrumental  = selfTaught ")

wt_ability_schoolInstrumental_selfTaught

Linear hypothesis test:
schoolInstrumental - selfTaught = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    490 2847.6                              
2    489 2815.1  1    32.483 5.6426 0.01791 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_schoolExtra_privateLesson <- linearHypothesis(lm_ability, "schoolExtracurricular  = privateLesson ")

wt_ability_schoolExtra_privateLesson

Linear hypothesis test:
schoolExtracurricular - privateLesson = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2830.3                           
2    489 2815.1  1     15.17 2.6351 0.1052

Code
wt_ability_schoolExtral_selfTaught <- linearHypothesis(lm_ability, "schoolExtracurricular  = selfTaught ")

wt_ability_schoolExtral_selfTaught

Linear hypothesis test:
schoolExtracurricular - selfTaught = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    490 2819.7                           
2    489 2815.1  1    4.6151 0.8017  0.371

Code
wt_ability_private_selfTaught <- linearHypothesis(lm_ability, "privateLesson  = selfTaught ")

wt_ability_private_selfTaught

Linear hypothesis test:
privateLesson - selfTaught = 0

Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender

  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    490 2865.0                                
2    489 2815.1  1    49.905 8.6688 0.003391 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald tests results indicated that the difference in coefficients between mindset of music ability and extracurricular school music club, mindset of music ability and private lesson, active engagement and extracurricular school music club, active engagement and self-teaching, school instrumental ensemble and self-teaching, and private lesson and self-teaching were not statistically significant, suggesting those predictors had the same effect on students’ global music self-concept. Coefficients for all other pairs were statistically different from each other.

Joint Significance (Wald Test)

Code
wt_ability_grade <- lmtest::waldtest(lm_ability, term = "grade")

wt_ability_grade
Wald test

Model 1: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + school + 
    gender
  Res.Df Df      F   Pr(>F)   
1    489                      
2    491 -2 5.4937 0.004371 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_school <- lmtest::waldtest(lm_ability, term = "school")

wt_ability_school
Wald test

Model 1: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    gender
  Res.Df Df      F   Pr(>F)   
1    489                      
2    491 -2 6.7475 0.001286 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_ability_gender <- lmtest::waldtest(lm_ability, term = "gender")

wt_ability_gender
Wald test

Model 1: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school + gender
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + 
    schoolExtracurricular + privateLesson + selfTaught + grade + 
    school
  Res.Df Df      F   Pr(>F)   
1    489                      
2    492 -3 3.9163 0.008797 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Scatter Plots

Code
#scatter plots with fitted regression lines
plot1_ability <- ggplot(data, aes(x = growth_mindset_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Mindset of Music Ability",
       y = "Ability Music Self-Concept") +
  theme_ipsum()

plot2_ability <- ggplot(data, aes(x = active_engagement_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Active Engagement",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot3_ability <- ggplot(data, aes(x = selfTaught_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Self-Teaching",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot4_ability <- ggplot(data, aes(x = schoolExtracurricular_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Extracurricular School Music Club",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot5_ability <- ggplot(data, aes(x = schoolInstrumental_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "School Instrumental Ensemble",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot6_ability <- ggplot(data, aes(x = privateLesson_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Private Lesson",
       y = "Global Music Self-Concept") +
  theme_ipsum()

#Combine plots into a single display (grid layout)
grid.arrange(plot1_ability, plot2_ability, plot3_ability, plot4_ability, plot5_ability, plot6_ability, ncol = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Alternative Scatter Plots Option

Code
#ALTERNATIVE TEST

#scatter plots with fitted regression lines
ggplot(data, aes(x = growth_mindset_z,  y = musical_ability_msc_z, color = school_music_elective)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Ability Music Self-Concept ~ Mindset of Music Ability",
       x = "Mindset of Music Ability",
       y = "Ability Music Self-Concept",
       color = "Enrolled in music elective?") +
  theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'

Code
ggplot(data, aes(x = active_engagement_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Active Engagement",
       y = "Global Music Self-Concept") +
  theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'

Code
ggplot(data, aes(x = selfTaught_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Self-Teaching",
       y = "Global Music Self-Concept") +
  theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'

Code
ggplot(data, aes(x = schoolExtracurricular_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Extracurricular School Music Club",
       y = "Global Music Self-Concept") +
  theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'

Code
plot5_ability <- ggplot(data, aes(x = schoolInstrumental_z,  y = musical_ability_msc_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "School Instrumental Ensemble",
       y = "Global Music Self-Concept") +
  theme_ipsum()

NON-ACADEMIC MSC - Stepwise Multiple Linear Regression

Intercept

Code
#Intercept
lm_msc_no_ability_step_0 <- lm(formula = msc_no_ability ~ 1, data = data)

lm_msc_no_ability_step_0

Call:
lm(formula = msc_no_ability ~ 1, data = data)

Coefficients:
(Intercept)  
      72.98  

Step 1

Code
#Step 1
lm_msc_no_ability_step_1 <- lm(formula = msc_no_ability ~ growth_mindset,  data = data)

summary_lm_msc_no_ability_step_1 <- summary(lm_msc_no_ability_step_1)

summary_lm_msc_no_ability_step_1

Call:
lm(formula = msc_no_ability ~ growth_mindset, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.961  -7.895   0.078   9.085  33.000 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     48.0790     2.5632  18.758   <2e-16 ***
growth_mindset   1.9869     0.1998   9.942   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.22 on 501 degrees of freedom
Multiple R-squared:  0.1648,    Adjusted R-squared:  0.1631 
F-statistic: 98.85 on 1 and 501 DF,  p-value: < 2.2e-16

Step 2

Code
#Step 2
lm_msc_no_ability_step_2 <- lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement,  data = data)

summary_lm_msc_no_ability_step_2 <- summary(lm_msc_no_ability_step_2)

summary_lm_msc_no_ability_step_2

Call:
lm(formula = msc_no_ability ~ growth_mindset + active_engagement, 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-28.6795  -5.5348   0.4917   5.7274  28.0049 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       22.28260    2.25993   9.860  < 2e-16 ***
growth_mindset     1.11784    0.15279   7.316 1.02e-12 ***
active_engagement  1.59886    0.07735  20.670  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.983 on 500 degrees of freedom
Multiple R-squared:  0.5496,    Adjusted R-squared:  0.5478 
F-statistic: 305.1 on 2 and 500 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_no_ability_1_2 <- anova(lm_msc_no_ability_step_1, lm_msc_no_ability_step_2)

f_test_result_lm_msc_no_ability_1_2
Analysis of Variance Table

Model 1: msc_no_ability ~ growth_mindset
Model 2: msc_no_ability ~ growth_mindset + active_engagement
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    501 74815                                  
2    500 40343  1     34472 427.24 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 3

Code
#Step 3
lm_msc_no_ability_step_3 <- lm(formula = msc_no_ability ~ growth_mindset  +
                   active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + 
                   outsideChoir + outsideInstrumental + 
                   privateLesson + selfTaught,  data = data)

summary_lm_msc_no_ability_step_3 <- summary(lm_msc_no_ability_step_3)

summary_lm_msc_no_ability_step_3

Call:
lm(formula = msc_no_ability ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.218  -5.703   0.358   5.803  28.814 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           23.716170   2.346261  10.108  < 2e-16 ***
growth_mindset         1.111693   0.164547   6.756 4.01e-11 ***
active_engagement      1.531543   0.083939  18.246  < 2e-16 ***
schoolChoir           -0.006283   0.236632  -0.027   0.9788    
schoolInstrumental    -0.427999   0.265330  -1.613   0.1074    
schoolExtracurricular -0.177074   0.306375  -0.578   0.5636    
outsideChoir           0.361198   0.453212   0.797   0.4259    
outsideInstrumental   -0.163464   0.515331  -0.317   0.7512    
privateLesson         -0.305774   0.220076  -1.389   0.1653    
selfTaught             0.764092   0.244716   3.122   0.0019 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.904 on 493 degrees of freedom
Multiple R-squared:  0.5636,    Adjusted R-squared:  0.5557 
F-statistic: 70.75 on 9 and 493 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_no_ability_2_3 <- anova(lm_msc_no_ability_step_2, lm_msc_no_ability_step_3)

f_test_result_lm_msc_no_ability_2_3
Analysis of Variance Table

Model 1: msc_no_ability ~ growth_mindset + active_engagement
Model 2: msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1    500 40343                              
2    493 39088  7    1254.8 2.2609 0.02845 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 4

Code
#Step 4
lm_msc_no_ability_step_4 <- lm(formula = msc_no_ability ~ growth_mindset  + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught  + grade + school + familyMusic + race + hispanic_latino + gender, data = data)

summary_lm_msc_no_ability_step_4 <- summary(lm.beta(lm_msc_no_ability_step_4))

summary_lm_msc_no_ability_step_4

Call:
lm(formula = msc_no_ability ~ growth_mindset + active_engagement + 
    schoolChoir + schoolInstrumental + schoolExtracurricular + 
    outsideChoir + outsideInstrumental + privateLesson + selfTaught + 
    grade + school + familyMusic + race + hispanic_latino + gender, 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-27.8726  -5.4640   0.3922   5.7536  28.7744 

Coefficients:
                                Estimate Standardized Std. Error t value
(Intercept)                   24.4065909           NA  3.4535028   7.067
growth_mindset                 1.1744169    0.2399478  0.1727753   6.797
active_engagement              1.4671971    0.5921259  0.0898091  16.337
schoolChoir                   -0.0398351   -0.0053960  0.2490325  -0.160
schoolInstrumental            -0.3717361   -0.0449088  0.2787228  -1.334
schoolExtracurricular         -0.1361955   -0.0152693  0.3121567  -0.436
outsideChoir                   0.3230472    0.0230407  0.4661907   0.693
outsideInstrumental           -0.1079797   -0.0068968  0.5314533  -0.203
privateLesson                 -0.3145132   -0.0470852  0.2293146  -1.372
selfTaught                     0.7664320    0.1119425  0.2497520   3.069
grade7                        -0.1127632   -0.0039222  1.0451684  -0.108
grade8                         0.1171442    0.0042469  1.0440654   0.112
schoolB                       -0.2212530   -0.0082851  1.0264244  -0.216
schoolC                       -0.8822040   -0.0292941  1.1601239  -0.760
familyMusicYes                -0.3262705   -0.0121840  0.8437192  -0.387
raceBlack or African American  1.9593949    0.0406179  2.6207337   0.748
raceSome Other Race            2.0222003    0.0296097  2.9907924   0.676
raceTwo or More Races          1.8872029    0.0364543  2.7166139   0.695
raceWhite                      1.1917803    0.0376185  2.2274437   0.535
hispanic_latinoYes             0.5822449    0.0093004  2.0399304   0.285
genderMale                    -1.8684377   -0.0698732  0.8824830  -2.117
genderNon-Binary              -0.0460634   -0.0003066  4.6226558  -0.010
genderPrefer not to say       -3.1897261   -0.0299035  3.2713764  -0.975
                              Pr(>|t|)    
(Intercept)                   5.60e-12 ***
growth_mindset                3.17e-11 ***
active_engagement              < 2e-16 ***
schoolChoir                    0.87298    
schoolInstrumental             0.18293    
schoolExtracurricular          0.66281    
outsideChoir                   0.48868    
outsideInstrumental            0.83908    
privateLesson                  0.17085    
selfTaught                     0.00227 ** 
grade7                         0.91413    
grade8                         0.91071    
schoolB                        0.82942    
schoolC                        0.44737    
familyMusicYes                 0.69915    
raceBlack or African American  0.45504    
raceSome Other Race            0.49928    
raceTwo or More Races          0.48759    
raceWhite                      0.59287    
hispanic_latinoYes             0.77544    
genderMale                     0.03475 *  
genderNon-Binary               0.99205    
genderPrefer not to say        0.33003    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.961 on 480 degrees of freedom
Multiple R-squared:  0.5697,    Adjusted R-squared:   0.55 
F-statistic: 28.89 on 22 and 480 DF,  p-value: < 2.2e-16
Code
# Calculate model fit improvement
f_test_result_lm_msc_no_ability_3_4 <- anova(lm_msc_no_ability_step_3, lm_msc_no_ability_step_4)

f_test_result_lm_msc_no_ability_3_4
Analysis of Variance Table

Model 1: msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught
Model 2: msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    493 39088                           
2    480 38543 13    544.64 0.5217 0.9116

Results Table

Code
#RESULTS TABLE Non-Academic MSC

# Extract coefficients and statistics
results_lm_msc_no_ability <- bind_rows(broom::tidy(lm.beta(lm_msc_no_ability_step_4), conf.int = TRUE)
)


# Round output to two decimal points
results_lm_msc_no_ability$estimate <- round(results_lm_msc_no_ability$estimate,2)

results_lm_msc_no_ability$std_estimate <- round(results_lm_msc_no_ability$std_estimate,2)

results_lm_msc_no_ability$std.error <- round(results_lm_msc_no_ability$std.error,2)

results_lm_msc_no_ability$statistic <- round(results_lm_msc_no_ability$statistic,2)

results_lm_msc_no_ability$p.value <- round(results_lm_msc_no_ability$p.value,3)

results_lm_msc_no_ability$conf.low <- round(results_lm_msc_no_ability$conf.low,2)

results_lm_msc_no_ability$conf.high <- round(results_lm_msc_no_ability$conf.high,2)


# Combine conf.low and conf.high into a single column called "95% CI" and change column names
results_lm_msc_no_ability <- results_lm_msc_no_ability %>%
  mutate("95% CI" = paste0("[", conf.low, ", ", conf.high, "]")) %>%
  select(-conf.low, -conf.high) %>%   # Remove the original conf.low and conf.high columns
  rename(b = estimate, 
         β = std_estimate,
           SE = std.error,
         t = statistic,
         p = p.value)

# Change order of columns
results_lm_msc_no_ability <- results_lm_msc_no_ability %>% 
  select(term, b, SE, β, t, "95% CI", p)


# Change output for p values undeer .001
results_lm_msc_no_ability <- results_lm_msc_no_ability %>% 
  mutate(p = ifelse(p < 0.001, "<.001", p))

# Change predictors' names
results_lm_msc_no_ability <- results_lm_msc_no_ability %>% 
  mutate(term = 
           case_match(term, 
                      "(Intercept)" ~ "(Intercept)",
                      "growth_mindset" ~ "Growth mindset", 
                      "active_engagement" ~ "Active Engagement",
                      "schoolChoir" ~ "School choir",
                      "schoolInstrumental" ~ "School instrumental",
                      "schoolExtracurricular" ~ "School music club",
                      "outsideChoir" ~ "Choir (outside school)",
                      "outsideInstrumental" ~ "instrumental (outside school)",
                      "privateLesson" ~ "Private lesson",
                      "selfTaught" ~ "Self-teaching",
                      "grade7" ~ "Grade 7",
                      "grade8" ~ "Grade 8",
                      "schoolB" ~ "School B",
                      "schoolC" ~ "School C",
                      "familyMusicYes" ~ "Family music engagement",
                      "raceBlack or African American" ~ "Black or African American",
                      "raceSome Other Race" ~ "Other race",
                      "raceTwo or More Races" ~ "Two or more races",
                      "raceWhite" ~ "White",
                      "hispanic_latinoYes" ~ "Ethnicity",
                      "genderMale" ~ "Male",
                      "genderNon-Binary" ~ "Non-Binary",
                      "genderPrefer not to say" ~ "Gender not disclosed")
                      )

print(results_lm_msc_no_ability, n= 100)
# A tibble: 23 × 7
   term                              b    SE     β     t `95% CI`      p    
   <chr>                         <dbl> <dbl> <dbl> <dbl> <chr>         <chr>
 1 (Intercept)                   24.4   3.45 NA     7.07 [NA, NA]      <.001
 2 Growth mindset                 1.17  0.17  0.24  6.8  [-0.1, 0.58]  <.001
 3 Active Engagement              1.47  0.09  0.59 16.3  [0.42, 0.77]  <.001
 4 School choir                  -0.04  0.25 -0.01 -0.16 [-0.49, 0.48] 0.873
 5 School instrumental           -0.37  0.28 -0.04 -1.33 [-0.59, 0.5]  0.183
 6 School music club             -0.14  0.31 -0.02 -0.44 [-0.63, 0.6]  0.663
 7 Choir (outside school)         0.32  0.47  0.02  0.69 [-0.89, 0.94] 0.489
 8 instrumental (outside school) -0.11  0.53 -0.01 -0.2  [-1.05, 1.04] 0.839
 9 Private lesson                -0.31  0.23 -0.05 -1.37 [-0.5, 0.4]   0.171
10 Self-teaching                  0.77  0.25  0.11  3.07 [-0.38, 0.6]  0.002
11 Grade 7                       -0.11  1.05  0    -0.11 [-2.06, 2.05] 0.914
12 Grade 8                        0.12  1.04  0     0.11 [-2.05, 2.06] 0.911
13 School B                      -0.22  1.03 -0.01 -0.22 [-2.03, 2.01] 0.829
14 School C                      -0.88  1.16 -0.03 -0.76 [-2.31, 2.25] 0.447
15 Family music engagement       -0.33  0.84 -0.01 -0.39 [-1.67, 1.65] 0.699
16 Black or African American      1.96  2.62  0.04  0.75 [-5.11, 5.19] 0.455
17 Other race                     2.02  2.99  0.03  0.68 [-5.85, 5.91] 0.499
18 Two or more races              1.89  2.72  0.04  0.69 [-5.3, 5.37]  0.488
19 White                          1.19  2.23  0.04  0.54 [-4.34, 4.41] 0.593
20 Ethnicity                      0.58  2.04  0.01  0.29 [-4, 4.02]    0.775
21 Male                          -1.87  0.88 -0.07 -2.12 [-1.8, 1.66]  0.035
22 Non-Binary                    -0.05  4.62  0    -0.01 [-9.08, 9.08] 0.992
23 Gender not disclosed          -3.19  3.27 -0.03 -0.98 [-6.46, 6.4]  0.33 

Overall Significance

Code
#overall significance

f_test_result_lm_msc_no_ability_0_4 <- anova(lm_msc_no_ability_step_0, lm_msc_no_ability_step_4)

print(f_test_result_lm_msc_no_ability_0_4) 
Analysis of Variance Table

Model 1: msc_no_ability ~ 1
Model 2: msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + 
    schoolInstrumental + schoolExtracurricular + outsideChoir + 
    outsideInstrumental + privateLesson + selfTaught + grade + 
    school + familyMusic + race + hispanic_latino + gender
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    502 89576                                  
2    480 38543 22     51032 28.888 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post Hoc Analyses

Equality of coefficients (Wald Test)

Code
#Visually compare coefficients
library (jtools)
library(broom.mixed)

lm <- lm(formula = msc_no_ability ~ growth_mindset + active_engagement + selfTaught + gender, data = data)

jtools::plot_summs(lm)


Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.

Code
library(car)

wt_mindset_active <- linearHypothesis(lm, "growth_mindset  = active_engagement ")

wt_mindset_active

Linear hypothesis test:
growth_mindset - active_engagement = 0

Model 1: restricted model
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F Pr(>F)  
1    497 39461                             
2    496 39157  1    303.43 3.8435 0.0505 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code
wt_mindset_selfTaught <- linearHypothesis(lm, "growth_mindset  = selfTaught ")

wt_mindset_selfTaught

Linear hypothesis test:
growth_mindset - selfTaught = 0

Model 1: restricted model
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1    497 39291                           
2    496 39157  1    133.37 1.6893 0.1943

Code
wt_active_selfTaught <- linearHypothesis(lm, "active_engagement  = selfTaught ")

wt_active_selfTaught

Linear hypothesis test:
active_engagement - selfTaught = 0

Model 1: restricted model
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught + 
    gender

  Res.Df   RSS Df Sum of Sq      F  Pr(>F)   
1    497 39843                               
2    496 39157  1    685.21 8.6795 0.00337 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald tests results indicated statistical difference between the coefficients of mindset of music ability and self-teaching, suggesting they had the same effect on students’ global music self-concept. Coefficients between the dyads mindset of music ability/active engagement and active engagement/self-teaching were statistically different from each other.

Joint Significance (Wald Test)

Code
wt_gender <- lmtest::waldtest(lm, term = "gender")

print(wt_gender)
Wald test

Model 1: msc_no_ability ~ growth_mindset + active_engagement + selfTaught + 
    gender
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught
  Res.Df Df      F Pr(>F)
1    496                 
2    499 -3 2.0287  0.109

Scatter Plots

Code
#scatter plots with fitted regression lines
plot1_no_ability <- ggplot(data, aes(x = growth_mindset_z,  y = msc_no_ability_z)) + 
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Mindset of Music Ability",
       y = "Non-Academic Music Self-Concept") +
  theme_ipsum()

plot2_no_ability <- ggplot(data, aes(x = active_engagement_z,  y = msc_no_ability_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Active Engagement",
       y = "Global Music Self-Concept") +
  theme_ipsum()

plot3_no_ability <- ggplot(data, aes(x = selfTaught_z,  y = msc_no_ability_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE)  +
  labs(x = "Self-Teaching",
       y = "Global Music Self-Concept") +
  theme_ipsum()

#Combine plots into a single display (grid layout)
grid.arrange(plot1_no_ability, plot2_no_ability, plot3_no_ability, ncol = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'